Subwavelength optical spatial solitons and three-dimensional localization in 
disordered ferroelectrics: towards metamaterials of nonlinear origin 
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We predict the existence of a novel class of multidimensional light localizations in out-of-equilibrium 
ferroelectric crystals. In two dimensions, the non-diffracting beams form at arbitrary low power level 
and propagate even when their width is well below the optical wavelength. In three dimensions, 
a novel form of subwavelength light bullets is found. The effects emerge when compositionally 
disordered crystals are brought to their metastable glassy state, and can have a profound impact 
on super-resolved imaging and ultra-dense optical storage, while resembling many features of the 
so-called metamaterials, as the suppression of evanescent waves. 



Introduction — Out-of-equilibrium materials display 
remarkable features, most of them still to be understood. 
Recent experiments in supercooled photorefractive crys- 
tals have allovifed the observation of "scale-free" optical 
solitons [ll supported by an extremely weak diffusive non- 
linearity [2|-y], which becomes active through the emer- 
gence of a dipolar glass with anomalously enhanced sus- 
ceptibility. These nonlinear beams have a truly remark- 
able feature: they are independent of size and intensity. 
Size-independence spurns a very basic and potentially 
ground-breaking exploration: are these scale-free optical 
beams capable of propagating even when their size is no- 
ticeably smaller than the optical wavelength? 

In this Letter we predict that glassy photorefractive 
ferroelectrics 0^ support a novel kind of multidimen- 
sional light localization at scales below the optical wave- 
length. The effect requires a huge susceptibility that can 
be harnessed in the out-of-equilibrium, or non-ergodic, 
phase, by acting on the previous history of the sample 
[l| . The finding is thus part of a still infant field of inves- 
tigation that focuses on using out-of-equilibrium optical 
materials to achieve novel and hereto unexplored effects 
0,111, a non-ergodic nonlinear optics that is rooted in sta- 
tistical mechanics, material science, and nonlinear wave 
propagation. Subwavelength propagation overcomes a 
basic limit to optical imaging and microscopy, i.e., that 
a light field can only propagate components of its spatial 
spectrum within the diffraction limit. High frequency 
components that correspond to details comparable to and 
smaller than the wavelength normally form evanescent 
waves that simply do not propagate. In stark contrast, 
in our predictions light leads to self-trapped beams of 
arbitrary intensity and widths, for which no diffraction 
limit holds. These predictions are based on an intensity- 
independent nonlinearity, which can be interpreted as the 
consequence of a widely tunable refractive index accom- 
panied by a transmission of evanescent fields. Put differ- 
ently, a nonlinearity-based meta-material 13" 14 
Intensity-independent PNR nonlinearity — Our model 



system is a compositionally disordered and impurity- 
doped photorefractive relaxor ferroelectric (e.g., the 
KLTN [5|). When it is rapidly cooled below the char- 
acteristic Bruns temperature Tb, it exhibits polar nano- 
regions (PNR). These are highly polarizable randomly- 
distributed ferroelectric-like regions that form a dipolar 
glass and provide an enormous enhancement of the pho- 
torefractive nonlinear optical response but retain very 
limited scattering losses f]\. The clue to the whole 
matter, originally discussed by Burns [ll], is that the 
optically-induced index of refraction change depends on 
the spatial average of the square of the crystal polariza- 
tion P. When the crystals display a zero polarization 
due to disorder (i.e., (P) ~ 0), averaging over disorder 
leads to a non-vanishing effect that depends exclusively 
on the mean-square of the polarization (i.e., (|Pp)). Op- 
tical response is thus directly correlated to the underlying 
nature of the crystal fluctuations, that become anoma- 
lously large and dependent on the crystal history in the 
out-of-equilibrium state. Specifically, averaging over the 
randomly oriented PNR [y, 0], the index of refraction 
perturbation is 
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where no is the bulk refractive index of the isotropic 
(disordered) crystals, g is the relevant component of 
the second-order electrooptic tensor, xpnr is the low- 
frequency electric response due to the PNR. The low- 
frequency electric field Edc is the space-charge field 
expressed in terms of the optical intensity / as P, 0| 
Edc = ~iKBT/q)\VI\/I, with T the crystal temper- 
ature, Kb the Boltzmann constant and q the elementary 
charge. Note that it is this dependence on / that bestows 
on all effects their characteristic intensity independence. 
Non-ergodicity implies that xpnr will depend on the his- 
tory of the sample, so that the same crystal at the same 
temperature will display radically different nonlinear op- 
tical responses [ij. Equation ([T]) holds because xpnr is 
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several orders of magnitudes greater than the susceptibil- 
ity of the crystals in the paraelectric phase xp, so that the 
nonlinear effect is mainly due to the PNR. Terms in the 
index perturbation that depend on the spontaneous po- 
larization Pq (i.e., in (Pg )) do not depend on the optical 
field, are negligible compared to the effects of the PNR, 
and are henceforth dropped. Finally, off-diagonal index 
modulation terms, that would amount to polarization ro- 
tation effects, are averaged out by the PNR disorder, and 
hence are negligible. 

Scale-free self-trapped paraxial beams — In the paraxial 
approximation, for a linearly polarized beam, the slowly- 
varying optical field A = / is the optical intensity) 
obeys the nonlinear equation 
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A = 0, (2) 



where we have introduced the characteristic length 
L = 4TTnleoy^XPNB.{KBT/q). 



(3) 



with g > (for commonly adopted ferroelectrics) and 
k = uino/c. The resulting model admits analytical self- 
trapped solutions, originally found in Q and experi- 
mentally investigated in [ij if condition L > A is ful- 
filled. For L = X one has the exact Gaussian solution 
A = aexp{—ij3z) with 



Aq exp 



(4) 



with /? = 2/kwQ, the nonlinear correction to the wave- 
vector k. Remarkably, in Eq.Q the waist wq of the soli- 
ton and its amplitude a-re free independent parame- 
ters. When L > X a, notable effect is that self-trapped 
solutions are given hy A = ae'^'^ with 
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and Ao and wo arbitrary constants (i.e., the "existence 
curve" is flat [2] ) , while 



(6) 



Interestingly, as L grows the beam loses its radial sym- 
metry, developing a square like profile. In Fig.dl]) we 
compare the two solutions. 

Suppression of evanescent waves and the absence of a 
diffraction limit — We now consider scale-free solutions 
for beam waists comparable and smaller than the wave- 
length. Specifically, we use the Helmholtz equation, 
which generalizes the paraxial equation ([2]) . The only ap- 
proximation we implement with respect to the Maxwell 
equations is neglecting coherent vectorial coupling (e.g.. 
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FIG. 1: Scale free solutions for various values of a = 
(L^/SX^). We show the profile in (a) of the Gaussian solution 
{a = 0.125) and in (b,c) of the generalized solution Eq.([5ll 



the transfer of energy from a linearly polarized beam to 
the orthogonal polarization): this is expected to play a 
negligible role for isotropic disordered crystals and for 
the diffusive nonlinearity, which only depends on the lo- 
cal intensity and not on the beam polarization 0, Q . The 
Helmholtz model reads as {n = no + An) 



Vi.+ (-) E^O, 



(7) 



whose propagation invariant solution is written as _B = 
aexjp{ikzz) where k^ is the overall wave- vector in the z 
direction (different from its nonlinear perturbation /3 in 
the paraxial model above). The scale- free Gaussian ^ 
is also solution of Eq. ([7|) with arbitrary amplitude Aq 
and waist wq, with 



(8) 



when the condition L = A is satisfied. Note that this 
solution exists (i.e. k^ is real) as long as 

A 



Wo > 



Trno 



(9) 



that is, for a beam waist that is (within moltiplicative 
constants) greater than the wavelength A/no. 

The key point is that more general solutions exist for 
L > X, analogous to what occurs in the paraxial case. 
These are still given by Eq.([5]) where, however, the wavec- 
tor is 



VV A ; wl 



(10) 



The corresponding lower limit for the waist is hence 

A 



Wo > 7- 



nno 



(11) 



The factor 7 plays a role similar to a Lorentz contraction 
term in special relativity, even if with a different meaning. 
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Specifically, for L > \/2A the lower limit for the waist is 
scaled by a factor 7 < 1, and correspondingly beams 
with a width smaller than the wavelength can propagate 
in the medium without distortion. Note that for this 
solution the evanescent waves are completely inhibited by 
the intensity independent scale-free nonlinearity, which 
leads to the conclusion that arbitrarily-low power beams 
with size below the wavelength can propagate. 
Non-paraxial regime — To investigate the formation 
of the two-dimensional (2D) soliton in the nonparaxial 
model we consider the forward propagating projection of 
the Helmholtz equation, which is written as 



(12) 



Equation (112]) reduces, under suitable limits, to the 
well-known unidirectional propagation equations (see 
IgI 17 [ and references therein) for the description of non- 
linear optics beyond the paraxial model. The basic dif- 
ference here is that the nonlinear refractive index An is 
retained under the square root, since our nonlinearity is 
intensity-independent and hence of the same order of the 
Laplacian even in the low intensity regions of the beam. 
After introducing the optical carrier with E = Aexjp{ikz) 
and the diffusive nonlinearity, the normalized dimension- 
less model reads as the non-paraxial normalized model 




1 + Jl + eV| - 2eo-(v • v) 



^ = 0, 



(13) 



where we introduce the dimensionless variables ^ — x/wq, 
77 — y/wo, and C = -z/zq with zq = kw^ the Rayleigh 
length. = A/Aq is the normalized optical field with 
Aq an arbitrary constant, and v = 'E,oc I {KbT / q) is the 
normalized space-charge field. In Eq.([2]) only two pa- 
rameters appear, the degree of non-par axiality e, which 
vanishes in the paraxial limit, and the strength of the 
scale- free nonlinearity ct (<t = 1/8 as L = A). 

Equation ([2]) can be numerically solved by Taylor ex- 
panding the square root, giving 



N 



= * E I [^L - 2a(v . v)] " V (14) 

n=l ^ ^ 

where N is the order of approximation selected (for a 
fixed e) in order to a have a given precision. 

To asses the existence of self-trapped scale-free soli- 
tons beyond the paraxial regime, we start by consider- 
ing the linear diffraction regime (cr = 0; L << A) in 
FiglSti., which shows the spreading of the waist of a Gaus- 
sian beam for an increasing order in the the solution of 
Eq. ()14p . For e > non-paraxial terms provide a more 
pronounced spreading if compared to the paraxial model 
(i.e., to = 1), and it is shown that for the case of 
e = 0.05 corrections become inconsequential for iV > 3. 
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FIG. 2: (a) Beam waist versus propagation for e = 0.05 for 
various orders of approximation of the nonparaxial equations 
(black line corresponds to standard paraxial models); as cr = 
the nonparaxial terms lead to a more pronounced spreading; 
for o = 0.125 an invariant propagation is attained at any non- 
paraxial order; (b) three dimensional comparison between the 
case (T = and a — 0.125 for = 3. 



In Figl^b we show the propagation of the beam for 
cr = 1/8 {L = A), including higher order diffraction. 
In such a nonlinear case, difFractionless propagation is 
achieved at any order and for any intensity. Note that 
the resulting beam is propagation invariant at any order 
N independently on the scale wq and on the intensity 
level, thus showing the fact that the proposed scale free 
solutions are indeed stable and exist for ultra-thin beams. 
Enhanced visibility — Next we numerically investigate 
conditions allowing super-resolution. We consider the 
propagation of two parallel beams for e = 0.01 (order 
N — 3) for various values of the ratio L/ X. For a given in- 
put pattern we find that there exist an optimal value cor- 
responding to mimimal input intensity distribution (i.e., 
image) distortion. In Figl^l we consider the propagation 
of two parallel beams (shown in FigI3^,c) for different 
values of a. We compare the (FiglS}^) linear propagation 
regime (cr = L/X = 0) and the case a — 0.15 (Figl^l;); 
for cr > (corresponding to high cooling rates see [ij) 
the visibility radically increases. In FiglSji we show the 
frange visibility versus cr, calculated as the ratio between 
the intensity peak value and the value at the center of 
the beam: the contrast increases (ideally diverges) when 
the visibility is higher. This result shows that even for 
beam sizes in the non-paraxial regime, at sufficiently high 
cooling rate, it is possible to propagate images with- 
out loss of information. We stress that these dynamics 
are attained at any intensity level, and hence represent 
a completely new regime in optical propagation, where 
the laws of diffraction and interference are largely modi- 
fied; the strength of nonlinearity being determined by a. 
For L >> A (cr >> 1/8), a Gaussian beam experiences 
a strong intensity-independent self-focusing, even below 
the non-paraxial regime (not reported). 
Three-dimensional subwavelength localization — A no- 
table property is the existence of three-dimensional lo- 
calized light bullets. Specifically the three-dimensional 
(3D) Hclhmholtz equation for the diffusive nonlinearity 
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7 A/no . The medium hence exhibits (for L >> A) an 
effective refractive index 



FIG. 3: Intensity independent super-resolution: (a,b) Prop- 
agation in the linear regime (L = cr = 0) of a double spot 
(panel a) in the non-paraxial model Ea. (|14|) at the third or- 
der {N = 3) with e = 0.01; (c,d) as in (a,b) with a = 0.2; (d) 
calculated visibility versus a. The initial image in panels a,c 
is the same and included for the sake of comparison. 



can be cast as a "nonlinear" eigenvalue problem: 



E 



(15) 



Eq. ([T5|) admits an exact 3D Gaussian solution A ~ 



Aq exp [- 



+ z^) /wq] when L = X, for any Aq 



and when wq — ^J'i/2T:\/nQ. In distinction to previous 
results, these solutions are not spatially scale-free and 
may only have a fixed waist, comparable with the wave- 
length in the medium. For L > A a solution exists (more 
general ones exist and will be reported elsewhere) and is 
given by 

A~Af^ cosh(\/2a;/wo) cos\Y{V2y / wo) cosh(-\/2z/wo) 

(16) 

with 7 as above (Eq.Q) and the waist wq — 
7-\/3/27rA/27rno. These 3D localized solutions have a 
waist smaller that the wavelength when 7 < 1 (L > 
■\/2A). They represent a novel form of light localization 
at any intensity level and with size comparable or smaller 
than the wavelength. We are not aware of other known 
light localizations, which can be described by an exact 
solution; it is a novel kind of bound state between the 
photo-induced charges and light which may be used to 
store information. 

Metamaterials of nonlinear origin — In standard optics, 
a Gaussian beam with waist wq has a spectral bandwidth 
of the order of I/wq, and the minimum waist such that 
the spectrum is contained in the Ewald circle (i.e., with- 
out evanescent waves) is given by A/no; in the scale- free 
regime here considered, due to the non-ergodic phase of 
glassy ferroelectrics, such a minimum waist is given by 



no 
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no^ >> no. (17) 



The PNR effect is therefore equivalent to an (inten- 
sity independent) refractive index, which can be largely 
tuned and increased, such that beams propagate without 
evanescent waves and without relevant scattering and ab- 
sorption losses. This shows that the specific nonlinearity 
here considered is able to provide those features that are 
the building blocks for the modern research on meta- 
materials, from a completely different perspective. 
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